{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\"> Registration Errors, Terminology and Interpretation <a href=\"https://mybinder.org/v2/gh/InsightSoftwareConsortium/SimpleITK-Notebooks/master?filepath=Python%2F68_Registration_Errors.ipynb\"><img style=\"float: right;\" src=\"https://mybinder.org/badge_logo.svg\"></a></h1>\n",
    "\n",
    "Registration is defined as the estimation of a geometric transformation aligning objects such that the distance between corresponding points on these objects is minimized, bringing the objects into alignment.\n",
    "\n",
    "When working with point-based registration algorithms we have three types of errors associated with our points (originally defined in [1]):  \n",
    "\n",
    "**Fiducial Localization Error (FLE)**: The error in determining the location of a point which is used to estimate the transformation. The most widely used FLE model is that of a zero mean Gaussian with independent, identically distributed errors. The figure below illustrates the various possible fiducial localization errors:\n",
    "\n",
    "<img src=\"fle.svg\" style=\"width:600px\"/><br><br>\n",
    "\n",
    "**Fiducial Registration Error (FRE)**: The error of the fiducial markers following registration, $\\|T(\\mathbf{p_f}) - \\mathbf{p_m}\\|$ where $T$ is the estimated transformation and the points $\\mathbf{p_f},\\;\\mathbf{p_m}$ were used to estimate $T$. \n",
    "\n",
    "**Target Registration Error (TRE)**: The error of the target fiducial markers following registration,$\\|T(\\mathbf{p_f}) - \\mathbf{p_m}\\|$ where $T$ is the estimated transformation and the points $\\mathbf{p_f},\\;\\mathbf{p_m}$ were **not** used to estimate $T$. \n",
    "\n",
    "\n",
    "Things to remember:\n",
    "1. TRE is the only quantity of interest, but in most cases we can only estimate its distribution.\n",
    "2. FRE should never be used as a surrogate for TRE as the TRE for a specific registration is uncorrelated with its FRE [2]. \n",
    "3. TRE is spatially varying.\n",
    "4. A good TRE is dependent on using a good fiducial configuration [3].\n",
    "5. The least squares solution to paired-point registration is sensitive to outliers.\n",
    "\n",
    "\n",
    "[1] \"[Registration of Head Volume Images Using Implantable Fiducial Markers](https://www.ncbi.nlm.nih.gov/pubmed/9263002)\", C. R. Maurer Jr. et al., *IEEE Trans Med Imaging*, 16(4):447-462, 1997.\n",
    "\n",
    "[2] \"[Fiducial registration error and target registration error are uncorrelated](http://spie.org/Publications/Proceedings/Paper/10.1117/12.813601)\", J. Michael Fitzpatrick, *SPIE Medical Imaging: Visualization, Image-Guided Procedures, and Modeling*, 7261:1–12, 2009.\n",
    "\n",
    "[3] \"[Fiducial point placement and the accuracy of point-based, rigid body registration](https://www.ncbi.nlm.nih.gov/pubmed/11322441)\", J. B. West et al., *Neurosurgery*, 48(4):810-816, 2001."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import SimpleITK as sitk\n",
    "import numpy as np\n",
    "import copy\n",
    "%matplotlib notebook\n",
    "from gui import PairedPointDataManipulation, display_errors\n",
    "import matplotlib.pyplot as plt\n",
    "from registration_utilities import registration_errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FLE, FRE, TRE empirical experimentation\n",
    "\n",
    "In the following cell you will use a user interface to experiment with the various concepts associated with FLE, FRE, and TRE. \n",
    "\n",
    "### Interacting with the GUI\n",
    "1. Change the mode and interact directly with the figure:\n",
    "    * Edit: Add pairs of fiducial markers with a left click and pairs of target markers with a right click.\n",
    "        * Markers can only be added prior to any manipulation (translation/rotation/noise/bias...)\n",
    "    * Translate: right-mouse button down + drag, anywhere in the figure.\n",
    "    * Rotate: right-mouse button down + drag, anywhere in the figure, will rotate the fiducials around their centroid  \n",
    "      (marked by a blue dot).\n",
    "2. Buttons:\n",
    "    * Clear Fiducials/Targets: Removes all fiducial/target marker pairs.\n",
    "    * Reset: Moves all marker pairs to the original locations.\n",
    "    * Noise: Add noise to the fiducial markers.\n",
    "    * Bias (FRE<TRE): Add bias to all fiducial markers so that FRE<TRE (move all markers in the same direction).\n",
    "    * Bias (FRE>TRE): Add bias to all fiducial markers so that FRE>TRE (move half of the markers in one direction and  \n",
    "      the other half in the opposite direction).\n",
    "    * Register: Align the two point sets using the paired fiducials. \n",
    "\n",
    "Marker glyphs:\n",
    "* Light red plus: moving fiducials\n",
    "* Dark red x: fixed fiducials\n",
    "* Light green circle: moving targets\n",
    "* Dark green square: fixed fiducials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "manipulation_interface = PairedPointDataManipulation(sitk.Euler2DTransform())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Sensitivity to outliers\n",
    "\n",
    "The least-squares solutions to the paired-point registration task used by the SimpleITK ``LandmarkBasedTransformInitializer`` method are optimal under the assumption of isotropic homogeneous zero mean Gaussian noise. They will readily fail in the presence of outliers (breakdown point is 0, a single outlier can cause arbitrarily large errors). \n",
    "\n",
    "The GUI above allows you to observe this qualitatively. In the following code cells we illustrate this quantitatively. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ideal_fixed_fiducials = [[23.768817532447077, 60.082971482049849], [29.736559467930949, 68.740980140058511],\n",
    "                   [37.639785274382561, 68.524529923608299], [41.994623984059984, 59.000720399798773]]\n",
    "ideal_fixed_targets = [[32.317204629221266, 60.732322131400501], [29.413978822769653, 56.403317802396167]]\n",
    "ideal_moving_fiducials = [[76.77857043206542, 30.557710579173616], [86.1401622129338, 25.76859196933914],\n",
    "                    [86.95501792478755, 17.904506579872375], [78.07960498849866, 12.346214284259808]]\n",
    "ideal_moving_targets = [[78.53588814928511, 22.166738486331596], [73.86559697098288, 24.481339720595585]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Registration with perfect data (no noise or outliers) \n",
    "fixed_fiducials = copy.deepcopy(ideal_fixed_fiducials)\n",
    "fixed_targets = copy.deepcopy(ideal_fixed_targets)\n",
    "moving_fiducials = copy.deepcopy(ideal_moving_fiducials)\n",
    "moving_targets = copy.deepcopy(ideal_moving_targets)\n",
    "\n",
    "# Flatten the point lists, SimpleITK expects a single list/tuple with coordinates (x1,y1,...xn,yn)\n",
    "fixed_fiducials_flat = [c for p in fixed_fiducials for c in p]\n",
    "moving_fiducials_flat = [c for p in moving_fiducials for c in p]\n",
    "\n",
    "transform = sitk.LandmarkBasedTransformInitializer(sitk.Euler2DTransform(), fixed_fiducials_flat, moving_fiducials_flat)\n",
    "\n",
    "FRE_information = registration_errors(transform, fixed_fiducials, moving_fiducials)\n",
    "TRE_information = registration_errors(transform, fixed_targets, moving_targets)\n",
    "FLE_values = [0.0]*len(moving_fiducials)\n",
    "FLE_information =  (np.mean(FLE_values), np.std(FLE_values), np.min(FLE_values), np.max(FLE_values), FLE_values) \n",
    "display_errors(fixed_fiducials, fixed_targets, FLE_information, FRE_information, TRE_information, title=\"Ideal Input\")\n",
    "\n",
    "# Change fourth fiducial to an outlier and register\n",
    "outlier_fiducial = [88.07960498849866, 22.34621428425981]\n",
    "FLE_values[3] = np.sqrt((outlier_fiducial[0] - moving_fiducials[3][0])**2 + \n",
    "                        (outlier_fiducial[1] - moving_fiducials[3][1])**2)\n",
    "moving_fiducials[3][0] = 88.07960498849866\n",
    "moving_fiducials[3][1] = 22.34621428425981\n",
    "\n",
    "moving_fiducials_flat = [c for p in moving_fiducials for c in p]\n",
    "\n",
    "transform = sitk.LandmarkBasedTransformInitializer(sitk.Euler2DTransform(), fixed_fiducials_flat, moving_fiducials_flat)\n",
    "\n",
    "FRE_information = registration_errors(transform, fixed_fiducials, moving_fiducials)\n",
    "TRE_information = registration_errors(transform, fixed_targets, moving_targets)\n",
    "FLE_information =  (np.mean(FLE_values), np.std(FLE_values), np.min(FLE_values), np.max(FLE_values), FLE_values) \n",
    "display_errors(fixed_fiducials, fixed_targets, FLE_information, FRE_information, TRE_information, title=\"Single Outlier\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FRE is not a surrogate for TRE\n",
    "\n",
    "In the next code cell we illustrate that FRE and TRE are not correlated. We first add the same fixed bias to \n",
    "all of the moving fiducials. This results in a large TRE, but the FRE remains zero. We then add a fixed bias to half of the moving fiducials and the negative of that bias to the other half. This results in a small TRE, but the FRE is large."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Registration with same bias added to all points\n",
    "fixed_fiducials = copy.deepcopy(ideal_fixed_fiducials)\n",
    "fixed_targets = copy.deepcopy(ideal_fixed_targets)\n",
    "moving_fiducials = copy.deepcopy(ideal_moving_fiducials)\n",
    "bias_vector = [4.5, 4.5]\n",
    "bias_fle = np.sqrt(bias_vector[0]**2 + bias_vector[1]**2)\n",
    "for fiducial in moving_fiducials:\n",
    "    fiducial[0] +=bias_vector[0]\n",
    "    fiducial[1] +=bias_vector[1]\n",
    "FLE_values = [bias_fle]*len(moving_fiducials)\n",
    "moving_targets = copy.deepcopy(ideal_moving_targets)\n",
    "\n",
    "# Flatten the point lists, SimpleITK expects a single list/tuple with coordinates (x1,y1,...xn,yn)\n",
    "fixed_fiducials_flat = [c for p in fixed_fiducials for c in p]\n",
    "moving_fiducials_flat = [c for p in moving_fiducials for c in p]\n",
    "\n",
    "transform = sitk.LandmarkBasedTransformInitializer(sitk.Euler2DTransform(), fixed_fiducials_flat, moving_fiducials_flat)\n",
    "\n",
    "FRE_information = registration_errors(transform, fixed_fiducials, moving_fiducials)\n",
    "TRE_information = registration_errors(transform, fixed_targets, moving_targets)\n",
    "FLE_information =  (np.mean(FLE_values), np.std(FLE_values), np.min(FLE_values), np.max(FLE_values), FLE_values) \n",
    "display_errors(fixed_fiducials, fixed_targets, FLE_information, FRE_information, TRE_information, title=\"FRE<TRE\")\n",
    "\n",
    "# Registration with bias in one direction for half the fiducials and in the opposite direction for the other half\n",
    "moving_fiducials = copy.deepcopy(ideal_moving_fiducials)\n",
    "pol = 1\n",
    "for fiducial in moving_fiducials:\n",
    "    fiducial[0] +=bias_vector[0]*pol\n",
    "    fiducial[1] +=bias_vector[1]*pol\n",
    "    pol*=-1.0\n",
    "FLE_values = [bias_fle]*len(moving_fiducials)\n",
    "moving_targets = copy.deepcopy(ideal_moving_targets)\n",
    "\n",
    "# Flatten the point lists, SimpleITK expects a single list/tuple with coordinates (x1,y1,...xn,yn)\n",
    "fixed_fiducials_flat = [c for p in fixed_fiducials for c in p]\n",
    "moving_fiducials_flat = [c for p in moving_fiducials for c in p]\n",
    "\n",
    "transform = sitk.LandmarkBasedTransformInitializer(sitk.Euler2DTransform(), fixed_fiducials_flat, moving_fiducials_flat)\n",
    "\n",
    "FRE_information = registration_errors(transform, fixed_fiducials, moving_fiducials)\n",
    "TRE_information = registration_errors(transform, fixed_targets, moving_targets)\n",
    "FLE_information =  (np.mean(FLE_values), np.std(FLE_values), np.min(FLE_values), np.max(FLE_values), FLE_values) \n",
    "display_errors(fixed_fiducials, fixed_targets, FLE_information, FRE_information, TRE_information, title=\"FRE>TRE\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fiducial Configuration\n",
    "\n",
    "Even when our model of the world is correct, no outliers and FLE is isotropic and homogeneous, the fiducial configuration\n",
    "has a significant effect on the TRE. Ideally you want the targets to be at the centroid of your fiducial configuration.\n",
    "\n",
    "This is illustrated in the code cell below. Translate, rotate and add noise to the fiducials, then register. The targets that are near the fiducials should have a better alignment than those far from the fiducials. \n",
    "\n",
    "Now, reset the setup. Where would you add two fiducials to improve the overall TRE? Experiment with various fiducial configurations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fiducials = [[31.026882048576109, 65.696247315510021], [34.252688500189009, 70.674602293864993], \n",
    "             [41.349462693737394, 71.756853376116084], [47.801075596963202, 68.510100129362826], \n",
    "             [52.47849495180192, 63.315294934557635]]\n",
    "targets = [[38.123656242124497, 64.397546016808718], [43.768817532447073, 63.748195367458059], \n",
    "           [26.833333661479333, 8.7698403891030861], [33.768817532447073, 8.120489739752438]]\n",
    "manipulation_interface = PairedPointDataManipulation(sitk.Euler2DTransform())\n",
    "manipulation_interface.set_fiducials(fiducials)\n",
    "manipulation_interface.set_targets(targets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FRE-TRE, and Occam's razor\n",
    "\n",
    "When we perform registration, our goal is to minimize the Target Registration Error. In practice it needs to be below a problem specific threshold for the registration to be useful.\n",
    "\n",
    "The target point(s) can be a single point or a region in space, and we want to minimize our registration error for this target. We go about this task by minimizing another quantity, in paired-point registration this is the FRE, in the case of intensity based registration we minimize an appropriate similarity metric. In both cases we expect that TRE is minimized indirectly.\n",
    "\n",
    "This can easily lead us astray, down the path of **overfitting**. In our 2D case, instead of using a rigid transformation with **three degrees of freedom** we may be tempted to use an affine transformation with **six degrees of freedom**. By introducing these additional degrees of freedom we will likely improve the FRE, but what about TRE?   \n",
    "\n",
    "In the cell below  you can qualitatively evaluate the effects of overfitting. Start by adding noise with no rotation or translation and then register. Switch to an affine transformation model and see how registration effects the fiducials and targets. You can then repeat this qualitative evaluation incorporating translation/rotation and noise. \n",
    "\n",
    "In this notebook we are working in an ideal setting, we know the appropriate transformation model is rigid. Unfortunately, this is often not the case. So which transformation model should you use? When presented with multiple competing hypotheses we select one using the principle of parsimony, often referred to as [Occam's razor](https://en.wikipedia.org/wiki/Occam%27s_razor). Our choice is to select the simplest model that can explain our observations. \n",
    "\n",
    "**In the case of registration, the transformation model with the least degrees of freedom that reduces the TRE below the problem specific threshold.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "simpleitk_error_allowed": "list index out of range"
   },
   "outputs": [],
   "source": [
    "fiducials = [[31.026882048576109, 65.696247315510021], \n",
    "             [41.349462693737394, 71.756853376116084], \n",
    "             [52.47849495180192, 63.315294934557635]]\n",
    "\n",
    "targets = [[38.123656242124497, 64.397546016808718], [43.768817532447073, 63.748195367458059]]\n",
    "manipulation_interface = PairedPointDataManipulation(sitk.Euler2DTransform())\n",
    "#manipulation_interface = PairedPointDataManipulation(sitk.AffineTransform(2))\n",
    "manipulation_interface.set_fiducials(fiducials)\n",
    "manipulation_interface.set_targets(targets)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
